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Abstract 

The Integrated Thermal Response Modeling System consists of a two- 
dimensional implicit thermal response and ablation program (TITAN), a commercialized 
finite-element thermal/mechanical analysis code (MARC), a high fidelity Navier-Stokes 
equation solver (GIANTS), and a boundary layer engineering correlation program 
(MEIT). This system is developed for fully coupled thermal response analysis of 
hypersonic entry vehicles with ablating heat-shield materials. The simulations performed 
by this integrated modeling system include flow' field, fluid/solid interaction, shape 
change, and heat-shield/structure material thermal response. Because of the severe 
aerothermal heating environment and the unique features of ablating heat-shield 
materials, none of the commercialized finite-element codes have the capability of 
simulating ablating heat-shield materials under hyper-velocity entry conditions. 

Therefore, TITAN was developed specifically for modeling ablating heat-shield. The 
structural materials underneath the heat-shield experience relatively benign heat loads, 
and are usually multi-dimensional and have complicated geometries. Thus, the 
commercialized finite-element code can be adapted for the thermal response simulation of 
structural materials. The approach taken in this work to perform full vehicle thermal 
analysis is simulating the ablating heat-shield using the TITAN code, and simulating the 
structural materials using MARC. The ablating heat-shield is treated as an outer boundary 
condition of structure. The continuity conditions of temperature and heat flux are 
imposed at the interface between TITAN and MARC. The aerothermal environments 
with fluid/solid interaction are predicted by coupling TITAN, GIANTS, and MEIT. The 
aerothermal environments and the thermal response solutions tor both ablating heat- 
shield and structure materials can be obtained simultaneously. Representative 
computations, including a proposed blunt body earth entry vehicle, a slender body re- 
entry vehicle, and an arc-jet test model are demonstrated. The results are presented and 
discussed in detail. 

Nomenclature 

B' = th i v i (\, , dimensionless mass blowing rate 

B, = pre-exponential constant in Eq.(S), s 1 

C, = skin friction coefficient 


C„ = Stanton number for heat transfer 

C M = Stanton number for mass transfer 

Cp = specific heat, J/kg-K 

E = internal energy, J/kg 

E, = activation energy in Eq.(8). J/kmol 

F = view factor 

h = enthalpy, J/kg 

H r = recovery enthalpy. J/kg 

k = thermal conductivity, W/m-K. 

m = mass flow rate, kg/nr-s 

P = pressure, N/nr 

q = heat flux, W/nr 

q ( = conductive heat flux. W/nr 

q R = radiative heat flux. W/nr 

Q = heat source term, J/m 

R = universal gas constant, J/kmol -K 

Rc = comer radius, cm 

Rn = nose radius, cm 

.v = stream length, m 

s = surface recession rate, m/s 

t = time, sec 

T = temperature, K 

u = velocity, m/s 

v 5 = species diffusion velocity, m/s 

v tt = surface nonrial velocity component, m/s 

Z* = coefficient in Eq.(3), defined in Ref. 14 

a = surface absorptance 

e = surface emissivity 

T = volume fraction of resin 

X = blowing reduction parameter 

0 = momentum boundary layer thickness, m 
p = total density, kg/nr 

a = Stefan-Boltzmann constant. W'nr-K. 4 
x = stress, N/ nr 
y = decomposition reaction order 
subscripts 
c = char 

e = boundary-layer edge 
g = pyrolysis gas 

1 = density component (A. B. and C) 
j = surface species 

m.n = index of computation cell 
v = virgin 


w = wall 


Introduction 

Spacecraft heat-shields often use Thermal Protection System (TPS) materials 
which pyrolyze and ablate at high temperature for mass-efficient rejection of acrothermal 
heat load. Pyrolysis is an internal decomposition of the solid w'hich releases gaseous 
species, whereas ablation is a combination of processes which consume heat-shield 
surface material. For design and sizing of ablating spacecraft TPS materials, it is 
imperative to have a reliable numerical procedure which can predict surface recession 
rate, in-depth pyrolysis, and internal temperature history. In addition to a numerical code, 
the thermal properties and the boundary' conditions of heat-shield must be precisely 
defined to perform an accurate simulation. The schematic diagram of Thermal Protection 
System is shown in Figure 1 . The front face of heat-shield is exposed to aerothermal 
heating generated by fluid/solid interactions. A high fidelity Computational Fluid 
Dynamics (CFD) code with appropriate surface thermal chemistry boundary conditions or 
an engineering code with semi-empirical correlations can be used to provide the 
aerothermal environments. The back face of heat-shield is attached to a structure. Thus, 
the structure must be included in the thermal response simulation to correctly predict the 
bond-line temperature history. The structure is usually with complicated geometries, but 
has relatively simple thermal response as compared with heat-shield. Hence, a 
commercialized finite-element material response code can easily be adapted for structural 
material thermal response simulation. 

In the past, a number of ablating heat-shield simulation programs w'ere developed. 
The Charring Material Thermal Response and Ablation Program 1 (CMA) was developed 
by the Aerotherm Corporation in the 1960’s. It solved the internal energy balance and 
decomposition equations coupled w'ith the ablating surface energy balance condition to 
simulate the response of ablative heat-shields in hypersonic flow's. The Fully Implicit 
Ablation and Thermal Response program 2 (FIAT) was developed by Chen and Milos at 
NASA Ames Research Center. The structure of FIAT is similar to that of One- 
dimensional Multi-Layer Implicit Thermal Response Simulation code, OMLITSv FIAT 
is numerically more stable, and solves much wdder range of problems compared with 
CMA, and has been used for TPS sizing calculations in various NASA's space missions. 
However, the governing equations solved in both CMA and FIAT codes are one- 
dimensional. Thus, these two codes are not applicable for conditions in which a one- 
dimensional assumption is not true, such as the nose-tip of a sharp body re-entry vehicle. 

The Ablation and Shape Change code 4 (ASCC) developed by the Aerotherm 
Corporation is an engineering code for fully coupled fluid/solid simulation with axi- 
symmetric or planar geometry. Nonetheless, the effect of pyrolysis gas due to internal 
decomposition w'as not implemented in ASCC, and its thermal diffusion equation w'as 
solved using an old fashioned finite-difference scheme. It has been shown that the 
accuracy of ASCC thermal response prediction is poor for both steady state and transient 
conditions'. The two-dimensional ablation code developed by Hogan et al at Sandia 
National Lab used the finite-control -volume method with unstructured grids". In this 
code, the mesh motion was modeled by assuming the mesh behaves as a linear elastic 



solid. In addition to two-dimensional thermal diffusion equation, a two-dimensional 
linear elastic solid mechanics model was applied to predict the motion of interior nodes. 
The effect of pyrolysis gas w'as not considered in Hogan’s w ork, and the interaction 
between solid and fluid was not discussed. 

In order to support the development of Ultra-High Temperature Ceramics 
(UHTCs) at NASA Ames Research Center, and to demonstrate their potential to enable 
sharp leading edges for hypersonic flight, an integrated analysis tool, MARC/ASCC\ was 
developed. Because of the unique structural, thermal, and chemical properties of UHTCs. 
the analysis requires full coupling of the flow' environment and the material thermal 
response. Thus, this integrated tool consists of a graphical user interface, a finite-element 
thermal and stress analysis package, a flow environment simulation code, and a material 
property database. The graphical user interface is know'n as Mentat, the finite-element 
material response code is MARC 7 , and the flow' environment simulation program is based 
on the Momentum and Energy Integral Technique Procedure. ME1T, in the ABRES 
Shape Change Code 4 . However, the integrated MARC/ASCC code is not capabe of 
simulating ablating TPS materials. 

Recently the Tw'o-dimensional Implicit Thermal Response and Ablation (TITAN) 
Program* w as developed to perform high fidelity thermal response and shape change 
simulation for charring materials. The governing equations include energy conservation 
and a three-component decomposition model w'ith a moving grid, and are discretized 
using a finite-volume approximation with general body fitted coordinates. The time 
accurate solution is achieved by the implicit time marching technique using Gauess- 
Siedel line relaxation with alternating sweeps 4 . 

A fully coupled fluid-solid simulation is required for an accurate prediction of 
shape change. This is because shape change can not be correctly determined without 
appropriate aerothermal heating, and aerothermal heating is very sensitive to the 
geometry of the solid surface. The Navier-Stokes solver 10 , GIANTS, and the Momentum 
and Energy Integral Technique Procedure 4 , MEIT, were successfully coupled with the 
TITAN code using a loosely coupled method for fluid/solid simulation. 

The purpose of this w r ork is to integrate the aerothermal environment code, 
ablating heat-shield simulation code, and structural material response code to perform full 
vehicle thermal response analysis. This system is built around the commercialized finite- 
element code. MARC. The TITAN, GIANTS, and MEIT codes are integrated with 
MARC through user supplied boundary condition subroutines. The ablating heat-shields 
are simulated using the TITAN code and the structural materials are simulated using 
MARC. This is because TITAN is a finite-volume code and thus is not as flexible as a 
finite-element code for solving complicated multi-dimensional objects. This integrated 
system can predict aerothermal heating, shape change, and heat-shield and structure 
thermal response simultaneously. To demonstrate the capability of this integrated 
material thermal response modeling system, the simulations of a proposed blunt body 
earth entry vehicle, a slender body re-entry vehicle, and an arc-jet test model arc 
performed. The results w ill be presented and discussed in the final paper. 
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Integrated Thermal Response Modeling System 

The flow chart for Integrated Thermal Response Modeling System is shown in 
Fig. 2. The four major components in this system are MARC. TITAN, GIANTS, and 
MEIT. The brief description of each code is listed in the following sections. The accuracy 
of codes has been previously studied 4 111 and is not the focus of this work. This modeling 
system is integrated around the MARC code. The front end of MARC is known as 
Mentat", which is a graphical user interface program. The material properties are 
provided from TPSX 12 , which is the TPS material property database developed at NASA 
Ames Research Center. The Multi-Component Ablation Thermochemistry code 1 ', MAT. 
is used to generate the tables of normalized mass blowing rate. 

MARC and TITAN are interfaced through the user provided heat flux boundary 
condition routine, FLUX. The continuity conditions of heat flux and temperature are 
imposed at the interface: 

Q MARC = ~ kVT 1 7/7 4 V 0 ) 

Tints =Imarc ( 2 ) 

At each time increment, the front-face heat flux boundary condition of MARC is updated 
using the temperature gradient calculated by TITAN (Eq. 1 ). and the back-face 
temperature boundary condition of TITAN is equal to the temperature computed by 
MARC (Eq. 2) at each interface point. 

TITAN and flow environment codes, including GIANTS and MEIT, are 
interfaced using the convective boundary condition. The conditions at the ablating surface 
are determined by convective and radiative heating and by surface thermochemical 
interactions with the boundary' layer gases. The surface energy balance equation 
employed is of the convective transfer coefficient type. This energy balance equation 
takes the following form 14 : 

Pc u ,£ h {H, W?r») + /? r w < C,j]T(Z* ( , )it'" - B'h»] + m,lh + mJn 

+ a„q R „-Fo£jt -q (K = 0 

The first term in Eq. 3 represents the sensible convective heat flux. The sum of the 
second, third, and fourth terms in Eq. 3 is defined as the total chemical energy at the 
surface. The Z* terms represent transport of chemical energy associated with chemical 
reactions at the wall and in the boundary layer 14 . The Z* driving forces for diffusive mass 
transfer include the effects of unequal diffusion coefficients. The fifth and sixth terms are 
the radiative heat fluxes absorbed and re-radiated by the wall, respectively, and the last 
term, q, w , represents the rate of conduction into material. Here B’ is the normalized mass 
blowing rate. The tables of B’ for charring materials can be generated using ACE 1 ' or 
MAT 1 ’" 

The computation starts from the MARC code. At each time increament, MARC 
calls TITAN to compute the new r heat flux boundary condition, and, at the same time 
TITAN obtains the back-face temperature boundary' condition from MARC and also calls 
the flow environment code to update the front-face convective boundary condition if 
necessary'. When the maximum local surface recession exceeds the predefined criterion 
since last surface convective heating was updated, a flow-field grid is generated based on 
the current body geometry, and then the flow simulation routine. GIANTS or MF.1T. is 



called to compute the aerothennal heating environment. Each call to the How 
environment routine is a steady stale calculation. In this work, the non-ablating cold wall 
heating is calculated by GIANTS using a 5-species air chemistry (N : , 0 : , NO, N, and O). 
and a blowing reduction parameter of 0.5 is used in TITAN to take into account the 
blockage due to mass blowing. As expected. GIANTS calculation is much more 
computationally intensive as than the thermal response computation. For a fully coupled 
high fidelity simulation, most of CPU time is consumed by the flow' environment 
calculation. Thus, the CPU time required for simulation is primarily determined by the 
efficiency of the GIANTS code. 

MARC 

The structure material thermal response and thermal stress analysis are performed 
using the commercialized finite-element code, MARC. Heat transfer and mechanical 
analysis can be done separately or fully coupled. The energy conservation law is 

P c/—-VkWT=Q (4) 

p rt 

The boundary conditions include prescribed temperature, heat flux, and convective heat 
transfer coefficient to the environment. For transient analysis, the initial temperature has 
to be specified to start the calculation. To perform thermal stress analysis, the MARC 
program uses an instantaneous thermal expansion coefficient defined as: 

de = a6 u dT (5) 

In most cases, the thermal expansion data is given with respect to a reference temperature 
T°, as: 

e ‘ h = a(T- T l ‘) (6) 

Hence, the conversion of the expansion data to the instantaneous thermal expansion 
coefficient becomes: 
da 

a = [a+ — (T-T k ')} (7) 

a I 

TITAN 

The ablating TPS material thermal response and shape change computation is 
performed using TITAN. The governing equations include energy conservation and three- 
component decomposition model. The surface energy balance condition (Eq. 3) is solved 
w'ith a moving grid to calculate the shape change due to surface recession. The internal 
energy balance is a transient thermal conduction equation with additional pyrolysis 
terms . 

“r 'T r 

pcr— = V-(kVT)-(h,-ii)V m v +m v -V/i, + spc r VT (8) 

c t 


The individual terms in Eq. 8 may be interpreted as follows: rate of storage of sensible 
energy, net rate of thermal conductive heal flux, pyrolysis energy-consumption rate, net 


rate of energy convected by pyrolysis, and convection rale of sensible energy due to 
coordinate system movement. 

A three-component decomposition model is used. The resin tiller is presumed to 
consist of two components which decompose separately, while the reinforcing material is 
the third component which can decompose. The instantaneous density of the composite is 
given by: 


p = \\p A +p H ) + {\-Y)p. (9) 

where A and B represent components of the resin, and C represents the reinforcing 
material. T is the volume fraction of resin and is an input quantity. Each of the three 
components can decompose following the relation: 


rt 


I. . C*m . ^ n , ' 1 1 / * , 'i 

= - B „, exp(— — )p„, ( ) + 5 Vp , 


RT 


P„ 


( 10 ) 


where p n is the residual or tenninal density of component i. and p,„ is the original density 
of component i ( i = A, B, and C). The motion of pyrolysis gas is assumed to be one- 
dimensional (along r| direction), and thus the mass flow rate of pyrolysis gas at the 
surface is calculated using the following approximation: 

1 ff u C P 

>*'=— f -(-i-)A„dn ( 11 ) 

/f h, dt 

Governing equations are discretized using a finite-volume approximation with a 
general body fitted coordinate system. The time accurate solution is achieved by the 
implicit time marching technique using Gauess-Siedel line relaxation with alternating 
sweeps. The computational grid is compressed during the course of computation to reflect 
the surface recession. 

MEIT 

The inviscid flow modeling in the MEIT code includes the prediction of the 
surface pressure distribution and the bow shock shape. The surface pressure ratio 
distribution is predicted by sophisticated semi-empirical correlations, each of w hich is 
specialized for specific regions of the body. The shock shape and shock standoff distance 
are predicted with an approximate integral solution of the continuity and momentum 
equations. 

The solution method for viscous boundary layer employed is known as the 
Momentum/Energy Integral Technique (MEIT). The following three basic equations are 
solved simultaneously: 

Integral Momentum Equation: 


1 <i , n 2 /lx ( > /’• « . HO dp 

_ ( /' p a () ) — 4 - 7 + ; 

/-"/■> u r ds ' 2 p ir p u ds 

Integral Energy Equation: 


(12) 


r"p, a, (//, , - /?„ ) ds 
Entrainment Relation: 


d h. - h f\ (h, , „ - ) 

-(/•>, »(A,, )</> ) = C h y f- + — — — — — 


(13) 


p a . »,.v N " =(2r)" /•// Re, -2" JrV M v> 


(14) 


Here, n = 1 for axi-symmetric fonnulation and n = 0 for 2-D planar flow. The inviscid 
flow solution provides the boundary layer edge conditions, and the surface energy 
balance computations provide the wall boundary conditions for the solution of the 
boundary layer integral equations. The boundary layer edge thermodynamic state is 
determined by lookup on pressure and entropy in a real gas Mollier table. Pressure is 
known from the inviscid flow solution, and entropy is calculated from consideration of 
the bow shock shape and boundary layer mass entrainment. The integral equations are 
solved in two ways: 1 .) series solutions near the stagnation point; 2.) finite difference 
numerical solutions away from the stagnation point. 

Various influence coefficients are defined in MEIT. The influence coefficients are 
factors derived by comparing the convective heat transfer to what it would be along the 
ideal plate for the same boundary layer state. The factors are derived in principle by 
considering one nonlinear behavioral mechanism at a time. Care is required to formulate 
the influence coefficients with the proper bases so that when the integral equations are 
integrated, the experimental data or exact solutions are accurately reproduced. The detail 
discussion regarding the implementation of influence coefficients are listed in Ref 4. 

The heating in the free molecular flow environment can be solved in closed 
analytical form. Prediction of the physics in the transitional regime presents a formidable 
challenge. The solutions of direct Monte-Carlo simulations require large amount of time. 
The approach taken in MEIT to compute transitional heating is very simple and consists 
of combing the free molecular and continuum value using an appropriate bridging 
relation. 


GIANTS 

The correlations implemented in MEIT are good for the slender body vehicle 
only. Thus, the flow simulation over a blunt body has to be performed using the GIANTS 
code. The GIANTS code solves the time-dependent conservation equations of mass, 
momentum, and energy for chemical and thermal non-equilibrium flowfield. The species 
mass conservation equation is given by 


+ — (A«,) = - 


■ ( P, V ,7 ) + V1 'v 


Op, 

Pt (\X , fix 

the momentum conservation is written as 


(15) 
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or 


(16) 
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and the energy conservation as 
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(17) 


The governing equations are discretized using finite-volume method. The numerical 
method that was used to solve these discretized equations is exactly the same as that used 
in the TITAN code. It's fully implicit and uses Gauss-Seidcl line relaxation. This 
technique has been shown to yield steady-state results efficiently. A 5-species air 
chemistry w'ith bifurcation diffusion model was implemented in this version of code. 


Results and Discussion 

Three test cases will be demonstrated in this work. The first case is a blunt body 
entry vehicle, the second one is a slender bode re-entry vehicle, and the third one is a fiat 
faced arc -jet model. Their temperature contours are shown in Figs 3 to 5. respectively. 
The detail of simulations will be discussed in the final paper. 
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